# older version of the orchaRd package
#devtools::install_github("itchyshin/orchard_plot", subdir = "orchaRd", force = TRUE, build_vignettes = TRUE)
#install.packages("cmdstanr", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))

pacman::p_load(tidyverse,
               purrr,
               metafor, 
               poolr,
               patchwork,
               orchaRd, # older version: 
               broom.mixed,
               here,
               nlme,
               pander,
               brms,
               kableExtra,
               formatR,
               knitr,
               brms,
               rstan,
               cmdstanr
               )

check_cmdstan_toolchain(fix = TRUE, quiet = TRUE)

# function to get what we need from these 2 models (you can include models in this function as well)
get_parmetersN<- function(i){
  
  # grand-mean centering of weights
  ln_c_weight <- scale(log(i[["weight"]]), center = TRUE, scale = TRUE)
  i[,"ln_c_weight"] <- ln_c_weight
  
  if(i[["nmeta"]][1] == 1 && i[["nstrain"]][1] == 1){
  
    # female model 
    model_f <- lme(log(data_point2) ~ sex*ln_c_weight, 
                   random = list(#metadata_group = ~ ln_c_weight, 
                                #strain_name = ~ ln_c_weight,
                                 date_of_experiment = ~ 1),
                   weights = varIdent(form = ~1 | sex),
                   control = lmeControl(opt = "optim"),
                   data = i)
    
    # male model
    model_m <- lme(log(data_point2) ~ relevel(sex, ref = "male")*ln_c_weight, 
                   random = list(#metadata_group = ~ ln_c_weight, 
                                 #strain_name = ~ ln_c_weight,
                                 date_of_experiment = ~ 1),
                   weights = varIdent(form = ~1 | sex),
                   control = lmeControl(opt = "optim"),
                   data = i)
    # neutral model
    model_n <- lme(log(data_point2) ~ sex*ln_c_weight, 
                   random = list(#metadata_group = ~ ln_c_weight, 
                                 #strain_name = ~ ln_c_weight,
                                 date_of_experiment = ~ 1),
                   #weights = varIdent(form = ~1 | sex),
                   control = lmeControl(opt = "optim"),
                   data = i)
      
  } else if (i[["nmeta"]][1] == 1) {
    
    # female model 
    model_f <- lme(log(data_point2) ~ sex*ln_c_weight, 
                   random = list(#metadata_group = ~ ln_c_weight, 
                     strain_name = ~ ln_c_weight,
                     date_of_experiment = ~ 1),
                   weights = varIdent(form = ~1 | sex),
                   control = lmeControl(opt = "optim"),
                   data = i)
    
    # male model
    model_m <- lme(log(data_point2) ~ relevel(sex, ref = "male")*ln_c_weight, 
                   random = list(#metadata_group = ~ ln_c_weight, 
                     strain_name = ~ ln_c_weight,
                     date_of_experiment = ~ 1),
                   weights = varIdent(form = ~1 | sex),
                   control = lmeControl(opt = "optim"),
                   data = i)
    # neutral model
    model_n <- lme(log(data_point2) ~ sex*ln_c_weight, 
                   random = list(#metadata_group = ~ ln_c_weight, 
                     strain_name = ~ ln_c_weight,
                     date_of_experiment = ~ 1),
                   #weights = varIdent(form = ~1 | sex),
                   control = lmeControl(opt = "optim"),
                   data = i)
    
  } else if (i[["nstrain"]][1] == 1){
    
    # female model 
    model_f <- lme(log(data_point2) ~ sex*ln_c_weight, 
                   random = list(metadata_group = ~ ln_c_weight, 
                     #strain_name = ~ ln_c_weight,
                     date_of_experiment = ~ 1),
                   weights = varIdent(form = ~1 | sex),
                   control = lmeControl(opt = "optim"),
                   data = i)
    
    # male model
    model_m <- lme(log(data_point2) ~ relevel(sex, ref = "male")*ln_c_weight, 
                   random = list(metadata_group = ~ ln_c_weight, 
                     #strain_name = ~ ln_c_weight,
                     date_of_experiment = ~ 1),
                   weights = varIdent(form = ~1 | sex),
                   control = lmeControl(opt = "optim"),
                   data = i)
    # neutral model
    model_n <- lme(log(data_point2) ~ sex*ln_c_weight, 
                   random = list(metadata_group = ~ ln_c_weight, 
                     #strain_name = ~ ln_c_weight,
                     date_of_experiment = ~ 1),
                   #weights = varIdent(form = ~1 | sex),
                   control = lmeControl(opt = "optim"),
                   data = i)
  } else {
    # female model 
    model_f <- lme(log(data_point2) ~ sex*ln_c_weight, 
                   random = list(metadata_group = ~ ln_c_weight, 
                                 strain_name = ~ ln_c_weight,
                                 date_of_experiment = ~ 1),
                   weights = varIdent(form = ~1 | sex),
                   control = lmeControl(opt = "optim"),
                   data = i)
    
    # male model
    model_m <- lme(log(data_point2) ~ relevel(sex, ref = "male")*ln_c_weight, 
                   random = list(metadata_group = ~ ln_c_weight, 
                                 strain_name = ~ ln_c_weight,
                                 date_of_experiment = ~ 1),
                   weights = varIdent(form = ~1 | sex),
                   control = lmeControl(opt = "optim"),
                   data = i)
    # neutral model
    model_n <- lme(log(data_point2) ~ sex*ln_c_weight, 
                   random = list(metadata_group = ~ ln_c_weight, 
                                 strain_name = ~ ln_c_weight,
                                 date_of_experiment = ~ 1),
                   #weights = varIdent(form = ~1 | sex),
                   control = lmeControl(opt = "optim"),
                   data = i)
  }
  # getting all we want
  females <- broom.mixed::tidy(model_f)
  males <- broom.mixed::tidy(model_m)
  # gets variance weights
  weights <- attr(model_f$modelStruct$varStruct, "weights")
  male_correction <- 1/weights[which(names(weights) == "male")[1]]
  female_correction <- 1/weights[which(names(weights) == "female")[1]]
  
  # get parameters
  parameter_name <- tolower(i[["parameter_name"]][1])
  procedure_name <- i[["procedure_name"]][1]# "procedure_name"
  m_n <- sum(i[["sex"]] == "male") # sample size for males 
  f_n <- sum(i[["sex"]] == "female") # N fo females
  f_intercept <- as.numeric(females[1, 4])
  f_intercept_se <- as.numeric(females[1, 5])
  f_slope <- as.numeric(females[3, 4])
  f_slope_se <- as.numeric(females[3, 5])
  m_intercept <- as.numeric(males[1, 4])
  m_intercept_se <- as.numeric(males[1, 5])
  m_slope  <- as.numeric(males[3, 4])
  m_slope_se  <- as.numeric(males[3, 5])
  fm_diff_int  <- as.numeric(males[2, 4])
  fm_diff_int_se  <- as.numeric(males[2, 5])
  fm_diff_int_p  <- as.numeric(males[2, 8])
  fm_diff_slope <- as.numeric(males[4, 4])
  fm_diff_slope_se <- as.numeric(males[4, 5])
  fm_diff_slope_p <- as.numeric(males[4, 8])
  
  # variance component
  #group_sd <- as.numeric(VarCorr(model_f)[,2][2])
  #g_slope_sd <- as.numeric(VarCorr(model_f)[,2][3])
  #batch_sd <- as.numeric(VarCorr(model_f)[,2][5])
  f_sd <- as.numeric(tail(VarCorr(model_f)[,2],1))*female_correction
  m_sd <- as.numeric(tail(VarCorr(model_f)[,2],1))*male_correction
  
  # model fit
  r_m <- sqrt(MuMIn::r.squaredGLMM(model_n)[1,1])
  r_c <- sqrt(MuMIn::r.squaredGLMM(model_n)[1,2])
  # putting it together
  paras <- data.frame(parameter_name, procedure_name, 
             f_n, m_n, f_intercept, f_intercept_se, f_slope, f_slope_se, 
             m_intercept, m_intercept_se, m_slope, m_slope_se, 
             fm_diff_int, fm_diff_int_se, fm_diff_int_p,
             fm_diff_slope, fm_diff_slope_se, fm_diff_slope_p,
             f_sd, m_sd, r_m, r_c)
  names(paras) <- c('parameter_name', 'procedure_name', 
                    'f_n', 'm_n','f_intercept', 'f_intercept_se', 'f_slope', 'f_slope_se',
                    'm_intercept', 'm_intercept_se', 'm_slope', 'm_slope_se',
                    'fm_diff_int', 'fm_diff_int_se', 'fm_diff_int_p',
                    'fm_diff_slope', 'fm_diff_slope_se', 'fm_diff_slope_p',
                    'f_sd', 'm_sd', 'r_m', 'r_c') # variance component
  invisible(paras)
  
}

# function to compare two models
get_strain_p<- function(i){
  
  # centering weights separately for each 
  ln_c_weight <- scale(log(i[["weight"]]), center = TRUE, scale = TRUE)
  i[,"ln_c_weight"] <- ln_c_weight
  
  if (i[["nmeta"]][1] == 1) {
    
    # model with strain as random factor 
    model_1 <- lme(log(data_point2) ~ sex*ln_c_weight, 
                   random = list(#metadata_group = ~ ln_c_weight, 
                     strain_name = ~ ln_c_weight,
                     date_of_experiment = ~ 1),
                   #weights = varIdent(form = ~1 | sex),
                   control = lmeControl(opt = "optim"),
                   data = i)
    
    # model without strain as random factor 
    model_2 <- lme(log(data_point2) ~ sex*ln_c_weight, 
                   random = list(#metadata_group = ~ ln_c_weight, 
                     #strain_name = ~ ln_c_weight,
                     date_of_experiment = ~ 1),
                   #weights = varIdent(form = ~1 | sex),
                   control = lmeControl(opt = "optim"),
                   data = i)
    
    
    
  }  else {

    # model with strain as random factor 
    model_1 <- lme(log(data_point2) ~ sex*ln_c_weight, 
                   random = list(metadata_group = ~ ln_c_weight, 
                                 strain_name = ~ ln_c_weight,
                                 date_of_experiment = ~ 1),
                   #weights = varIdent(form = ~1 | sex),
                   control = lmeControl(opt = "optim"),
                   data = i)
    # model without strain as random factor 
    model_2 <- lme(log(data_point2) ~ sex*ln_c_weight, 
                   random = list(metadata_group = ~ ln_c_weight, 
                                 #strain_name = ~ ln_c_weight,
                                 date_of_experiment = ~ 1),
                   #weights = varIdent(form = ~1 | sex),
                   control = lmeControl(opt = "optim"),
                   data = i)
    
  }
  
  # anova 
  p_value <- anova(model_1, model_2)$p[[2]]
  delta_aic <- anova(model_1, model_2)$AIC[1] - anova(model_1, model_2)$AIC[2] 
  
  # get parameters
  parameter_name <- tolower(i[["parameter_name"]][1])
  procedure_name <- i[["procedure_name"]][1]# "procedure_name"
  
  paras <- data.frame(parameter_name, procedure_name, delta_aic, p_value)
  names(paras) <- c('parameter_name', 'procedure_name', 'delta_aic', 'p_value') # variance component
  invisible(paras)
  
}


# getting ride of traits which do not run
get_para_poss <- possibly(.f = get_parmetersN, 
                          otherwise = NULL)


# getting ride of traits with do not run

get_para_poss2 <- possibly(.f = get_strain_p, 
                          otherwise = NULL)

# functions

orchard_plot2 <- function (object, mod = "Int", xlab, N = "none", alpha = 0.5, 
          angle = 90, cb = FALSE, k = TRUE, transfm = c("none", "tanh"), 
          point.size = 2.5, branch.size = 5,
          condition.lab = "Condition", legend.on = TRUE) 
{
  transfm <- match.arg(transfm)
  if (any(class(object) %in% c("rma.mv", "rma"))) {
    if (mod != "Int") {
      object <- mod_results(object, mod)
    }
    else {
      object <- mod_results(object, mod = "Int")
    }
  }
  mod_table <- object$mod_table
  data <- object$data
  data$moderator <- factor(data$moderator, levels = mod_table$name, 
                           labels = mod_table$name)
  data$scale <- (1/sqrt(data[, "vi"]))
  legend <- "Precision (1/SE)"
 
  # sample size 
  if(any(N != "none")){
    data$scale <- N
    legend <- "Sample size (N)" # we want to use italic
  }
  
  if (transfm == "tanh") {
    cols <- sapply(mod_table, is.numeric)
    mod_table[, cols] <- Zr_to_r(mod_table[, cols])
    data$yi <- Zr_to_r(data$yi)
    label <- xlab
  }
  else {
    label <- xlab
  }
  mod_table$K <- as.vector(by(data, data[, "moderator"], function(x) length(x[, 
                                                                              "yi"])))
  group_no <- length(unique(mod_table[, "name"]))
  cbpl <- c("#E69F00", "#009E73", "#F0E442", "#0072B2", "#D55E00", 
            "#CC79A7", "#56B4E9", "#999999")
  if (names(mod_table)[2] == "condition") {
    condition_no <- length(unique(mod_table[, "condition"]))
    plot <- ggplot2::ggplot() + ggbeeswarm::geom_quasirandom(data = data, 
                                                             ggplot2::aes(y = yi, x = moderator, size = scale, 
                                                                          colour = moderator), alpha = alpha) + ggplot2::geom_hline(yintercept = 0, 
                                                                                                                                    linetype = 2, colour = "black", alpha = alpha) + 
      ggplot2::geom_linerange(data = mod_table, ggplot2::aes(x = name, 
                                                             ymin = lowerCL, ymax = upperCL), size = branch.size, 
                              position = ggplot2::position_dodge2(width = 0.3)) + 
      ggplot2::geom_pointrange(data = mod_table, ggplot2::aes(y = estimate, 
                                                              x = name, ymin = lowerPR, ymax = upperPR, shape = as.factor(condition), 
                                                              fill = name), size = 0.5, position = ggplot2::position_dodge2(width = 0.3)) + 
      ggplot2::scale_shape_manual(values = 20 + (1:condition_no)) + 
      ggplot2::coord_flip() + ggplot2::theme_bw() + ggplot2::guides(fill = "none", 
                                                                    colour = "none") + ggplot2::theme(legend.position = c(0, 
                                                                                                                          1), legend.justification = c(0, 1)) + ggplot2::theme(legend.title = ggplot2::element_text(size = 9)) + 
      ggplot2::theme(legend.direction = "horizontal") + 
      ggplot2::theme(legend.background = ggplot2::element_blank()) + 
      ggplot2::labs(y = label, x = "", size = legend) + 
      ggplot2::scale_size_continuous(breaks = c(200, 2000, 20000), guide = guide_legend()) + 
      ggplot2::labs(shape = condition.lab) + ggplot2::theme(axis.text.y = ggplot2::element_text(size = 10, 
                                                                                                colour = "black", hjust = 0.5, angle = angle))
    plot <- plot + ggplot2::annotate("text", y = (max(data$yi) + 
                                                    (max(data$yi) * 0.1)), x = (seq(1, group_no, 1) + 
                                                                                  0.3), label = paste("italic(k)==", mod_table$K[1:group_no]), 
                                     parse = TRUE, hjust = "right", size = 3.5)
  }
  else {
    plot <- ggplot2::ggplot(data = mod_table, ggplot2::aes(x = estimate, 
                                                           y = name)) + ggbeeswarm::geom_quasirandom(data = data, 
                                                                                                     ggplot2::aes(x = yi, y = moderator, size = scale, 
                                                                                                                  colour = moderator), groupOnX = FALSE, alpha = alpha) + 
      ggplot2::geom_errorbarh(ggplot2::aes(xmin = lowerPR, 
                                           xmax = upperPR), height = 0, show.legend = FALSE, 
                              size = 0.5, alpha = 0.6) + ggplot2::geom_errorbarh(ggplot2::aes(xmin = lowerCL, 
                                                                                              xmax = upperCL), height = 0, show.legend = FALSE, 
                                                                                 size = branch.size) + ggplot2::geom_vline(xintercept = 0, 
                                                                                                                   linetype = 2, colour = "black", alpha = alpha) + 
      ggplot2::geom_point(ggplot2::aes(fill = name), size = point.size, 
                          shape = 21) + ggplot2::theme_bw() + ggplot2::guides(fill = "none", 
                                                                              colour = "none") + ggplot2::theme(legend.position = c(1, 
                                                                                                                                    0), legend.justification = c(1, 0)) + ggplot2::theme(legend.title = ggplot2::element_text(size = 9)) + 
      ggplot2::theme(legend.direction = "horizontal") + 
      ggplot2::theme(legend.background = ggplot2::element_blank()) + 
      ggplot2::labs(x = label, y = "", size = legend) + 
      ggplot2::scale_size_continuous(breaks = c(200, 2000, 20000), guide = guide_legend()) + 
      ggplot2::theme(axis.text.y = ggplot2::element_text(size = 10, 
                                                         colour = "black", hjust = 0.5, angle = angle))
    if (k == TRUE) {
      plot <- plot + ggplot2::annotate("text", x = (max(data$yi) + 
                                                      (max(data$yi) * 0.1)), y = (seq(1, group_no, 
                                                                                      1) + 0.3), label = paste("italic(k)==", mod_table$K), 
                                       parse = TRUE, hjust = "right", size = 3.5)
    }
  }
  if (cb == TRUE) {
    plot <- plot + ggplot2::scale_fill_manual(values = cbpl) + 
      ggplot2::scale_colour_manual(values = cbpl)
  }
  
  if (legend.on == FALSE){
    plot <- plot + ggplot2::theme(legend.position = "none")
  }
  
  
  return(plot)
}

# mod_result old

#' @title get_est
#' @description Function gets estimates from rma objects (metafor)
#' @param model rma.mv object
#' @param mod the name of a moderator. If meta-analysis (i.e. no moderator, se mod = "Int")
#' @author Shinichi Nakagawa - s.nakagawa@unsw.edu.au
#' @author Daniel Noble - daniel.noble@anu.edu.au
#' @export

get_est <- function (model, mod) {
  name <- firstup(as.character(stringr::str_replace(row.names(model$beta), {{mod}}, "")))
  
  estimate <- as.numeric(model$beta)
  lowerCL <- model$ci.lb
  upperCL <- model$ci.ub
  
  table <- tibble::tibble(name = factor(name, levels = name, labels = name), estimate = estimate, lowerCL = lowerCL, upperCL = upperCL)
  
  return(table)
}


#' @title get_pred
#' @description Function to get prediction intervals (crediblity intervals) from rma objects (metafor)
#' @param model rma.mv object
#' @param mod the name of a moderator
#' @author Shinichi Nakagawa - s.nakagawa@unsw.edu.au
#' @author Daniel Noble - daniel.noble@anu.edu.au
#' @export

get_pred <- function (model, mod) {
  
  name <- firstup(as.character(stringr::str_replace(row.names(model$beta), {{mod}}, "")))
  len <- length(name)
  
  if(len != 1){
    newdata <- matrix(NA, ncol = len, nrow = len)
    
    pred <- metafor::predict.rma(model, newmods = diag(len),
                                 tau2.levels = 1:len,
                                 gamma2.levels = 1:len)
  }
  else {
    pred <- metafor::predict.rma(model)
  }
  lowerPR <- pred$cr.lb
  upperPR <- pred$cr.ub
  
  table <- tibble::tibble(name = factor(name, levels = name, labels = name), lowerPR = lowerPR, upperPR = upperPR)
  return(table)
}

#' @title firstup
#' @description Uppercase moderator names
#' @param x a character string
#' @author Shinichi Nakagawa - s.nakagawa@unsw.edu.au
#' @author Daniel Noble - daniel.noble@anu.edu.au
#' @return Returns a character string with all combinations of the moderator level names with upper case first letters
#' @export
firstup <- function(x) {
  substr(x, 1, 1) <- toupper(substr(x, 1, 1))
  x
}

#' @title get_data
#' @description Collects and builds the data used to fit the rma.mv or rma model in metafor
#' @param model rma.mv object
#' @param mod the moderator variable
#' @author Shinichi Nakagawa - s.nakagawa@unsw.edu.au
#' @author Daniel Noble - daniel.noble@anu.edu.au
#' @return Returns a data frame
#' @export
#'
get_data <- function(model, mod){
  X <- as.data.frame(model$X)
  names <- vapply(stringr::str_split(colnames(X), {{mod}}), function(x) paste(unique(x), collapse = ""), character(1L))
  
  moderator <- matrix(ncol = 1, nrow = dim(X)[1])
  
  for(i in 1:ncol(X)){
    moderator <- ifelse(X[,i] == 1, names[i], moderator)
  }
  moderator <- firstup(moderator)
  yi <- model$yi
  vi <- model$vi
  type <- attr(model$yi, "measure")
  
  data <- data.frame(yi, vi, moderator, type)
  return(data)
  
}

#' @title mod_results
#' @description Using a metafor model object of class rma or rma.mv it creates a table of model results containing the mean effect size estimates for all levels of a given categorical moderator, their corresponding confidence intervals and prediction intervals
#' @param model rma.mv object
#' @param mod the name of a moderator; put "Int" if the intercept model (meta-analysis) or no moderators.
#' @return A data frame containing all the model results including mean effect size estimate, confidence and prediction intervals
#' @author Shinichi Nakagawa - s.nakagawa@unsw.edu.au
#' @author Daniel Noble - daniel.noble@anu.edu.au
#' @examples
#' \dontrun{data(eklof)
#' eklof<-metafor::escalc(measure="ROM", n1i=N_control, sd1i=SD_control,
#' m1i=mean_control, n2i=N_treatment, sd2i=SD_treatment, m2i=mean_treatment,
#' data=eklof)
#' # Add the unit level predictor
#' eklof$Datapoint<-as.factor(seq(1, dim(eklof)[1], 1))
#' # fit a MLMR - accouting for some non-independence
#' eklof_MR<-metafor::rma.mv(yi=yi, V=vi, mods=~ Grazer.type-1, random=list(~1|ExptID,
#' ~1|Datapoint), data=eklof)
#' results <- mod_results(eklof_MR, mod = "Grazer.type")
#' }
#' @export

mod_results <- function(model, mod) {
  
  if(all(class(model) %in% c("rma.mv", "rma.uni", "rma")) == FALSE) {stop("Sorry, you need to fit a metafor model of class rma.mv or rma")}
  
  data <- get_data(model, mod)
  
  # Get confidence intervals
  CI <- get_est(model, mod)
  
  # Get prediction intervals
  PI <- get_pred(model, mod)
  
  model_results <- list(mod_table = cbind(CI, PI[,-1]), data = data)
  
  class(model_results) <- "orchard"
  
  return(model_results)
  
}
# TODO - I think we can improve `mod` bit?

#' @title print.orchard
#' @description Print method for class 'orchard'
#' @param object x an R object of class orchard
#' @param ... Other arguments passed to print
#' @author Shinichi Nakagawa - s.nakagawa@unsw.edu.au
#' @author Daniel Noble - daniel.noble@anu.edu.au
#' @return Returns a data frame
#' @export
#'
print.orchard <- function(object, ...){
  return(object$mod_table)
}

# loading data
allometry <- readRDS(here("data/allometryNEW.rds"))


#STEP 1 remove rows with missing  data and NA 
allometrynew<-allometry[complete.cases(allometry),]

# getting rid of NA for data_point and weight

allometrynew2 <- allometrynew %>% 
  filter(!is.na(data_point), !is.na(weight)) %>% 
  group_by(parameter_name, sex, metadata_group, strain_name) %>%
  mutate(count = n()) %>% 
  ungroup() %>% 
  group_by(parameter_name) %>% # adjusting interval data
  mutate(min_val = min(data_point),
         data_point2 = if_else(min_val > 0, data_point, data_point + abs(min_val)),
         min_val2 = min( data_point[data_point!=min(data_point)]),
         data_point2 = if_else(min_val == 0, data_point2 + min_val2, data_point2),
         ratio_int =  if_else(min_val > 0, "ratio", "interval"),
         new_min = min(data_point2),
         nmeta = n_distinct(metadata_group),
         nstrain = n_distinct(strain_name),
         sex = as.factor(sex),
         parameter_name = if_else(parameter_name == "Latency to fall_Mean", 
                                  "Latency to fall mean"  , parameter_name)) %>% 
  ungroup() %>% 
  filter(count > 49) %>% # this can be adjusted
  filter(parameter_name != "BMC/Body weight", 
         parameter_name != "Body weight",  
         parameter_name != "Body Weight", 
         parameter_name != "Body weight after experiment" , 
         parameter_name != "Body weight before experiment",
         parameter_name != "Test duration") %>% 
  filter(!is.infinite(data_point2), !is.infinite(log(data_point2))) # removing infite and 0

# dim(allometry)
# dim(allometrynew)
# dim(allometrynew2)
# 
# # the number of traits
# length(unique(allometrynew2$parameter_name))
# 
# # the number of substrains
# length(unique(allometrynew2$strain_name))
# strain information
unique(allometrynew2$strain_name)

# # check there is no 0
# sum(is.infinite(log(allometrynew2$data_point2)))

# # the number of interval scale traits
# allometrynew2 %>% group_by(parameter_name) %>% summarise(ratio_int = ratio_int[1]) -> sum_ri
# sum(sum_ri$ratio_int == "interval")


#split dataframe by parameter to generate a list of dfs
#all_list<-split(allometrynew2, allometrynew2$parameter_name)

#saveRDS(all_list, file = here("data", "dat_list2.rds"))

# loading data
dat_list <- readRDS(here("data/dat_list2.rds"))

# grouping for category and parameter_group (this is from Zajitschek et al.
# 2020 eLife; slightly modfied)
dat_category <- read_csv(here("data/cateogry_parameter3.csv"))

#run function across list of matrices
processing <-map_dfr(dat_list, get_para_poss)
dat <- data.frame(processing, row.names = NULL)

dat %>% left_join(dat_category, by = ("parameter_name" = "parameter_name") ) %>% arrange(Category)  -> dat

dim(dat)

#write_csv(dat, here("data/test4.csv"))
write_csv(dat, here("data/data_parameters7.csv"))


# first getting p values - the contrasts between males and females for 

dat <-read_csv(here("data/data_parameters7.csv"))

#assess number of traits with sig shifts in intercept and slope

# getting lnVR to compare SDs and SD and Zr (variance for Zr)

dat %>% mutate(lnVR = log(f_sd/m_sd) + 1/(2*(f_n-3)) - 1/(2*(m_n-3)), 
               VlnVR = 1/(2*(f_n-3)) + 1/(2*(m_n-3)), 
               low_lnVR = lnVR - qnorm(0.975)*VlnVR, 
               high_lnVR = lnVR + qnorm(0.975)*VlnVR,
               t_val_sd = lnVR/sqrt(VlnVR),
               p_val_sd = 2*(1-pt(abs(t_val_sd), f_n-1 + m_n-1)),
               # r squared
               Zr = atanh(r_m),
               VZr = 1/((f_n + m_n) - 3)
               ) -> dat

write_csv(dat, here("data/data_parameters8.csv"))

# loading data
dat <- read_csv(here("data/data_parameters8.csv"))
# creating observation level random effect
dat$obs <- 1:dim(dat)[1]
# making character strings into factors
dat <- dat %>%
    mutate_if(is.character, as.factor)
# visualizing
kable(dat, "html") %>%
    kable_styling("striped", position = "left") %>%
    scroll_box(width = "100%", height = "250px")